# Replication Code for "Shocked into Service: Free Trade and the American South’s Military Burden," International Interactions
# Adam Dean & Jonathan Obert

# Import Data
y3 <- read.csv("ShockedIntoService.csv", header = TRUE)

# Table 3: Main Results

# Model 1

lm1.1 <- lm(Enrollment ~ lag1Enrollment
            + TradeLayoffs
            + Vetspc + Base + Rep2PartyVoteShare
            + EVANRATE + netMLS
            + median_age
            + lPopulation
            + IncomePC
            + Unemployment
            + Percent_Black
            + as.factor(Year)
            + State
            , data = 
              y3)
summary(lm1.1)

# Standard Errors clustered at county-level

cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y3,lm1.1,y3$GeoFips)

cl3.1 <- cl[, 1:4]
cl3.1 # Table 3, Model 1 results

# Model 2

lm1.2 <- lm(Enrollment ~ lag1Enrollment
            + TradeLayoffs*South
            + Vetspc + Base + Rep2PartyVoteShare
            + EVANRATE + netMLS
            + median_age
            + lPopulation
            + IncomePC
            + Unemployment
            + Percent_Black
            + as.factor(Year)
            #+ State
            , data = 
              y3)
summary(lm1.2)

# Standard Errors clustered at county-level

cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y3,lm1.2,y3$GeoFips)

cl3.2 <- cl[, 1:4]
cl3.2 # Table 3, Model 2 results

# Model 3

lm1.3 <- lm(Enrollment ~ lag1Enrollment
            + log(TradeLayoffs+1)*South
            + Vetspc + Base + Rep2PartyVoteShare
            + EVANRATE + netMLS
            + median_age
            + lPopulation
            + IncomePC
            + Unemployment
            + Percent_Black
            + as.factor(Year)
            #+ State
            , data = 
              y3)
summary(lm1.3)

# Standard Errors clustered at county-level
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y3,lm1.3,y3$GeoFips)

cl3.3 <- cl[, 1:4]
cl3.3 # Table 3, Model 3 results

# Table 4: What Makes the South Unique?

# 1. Veterans
lm1.1 <- lm(Enrollment ~ lag1Enrollment
            + TradeLayoffs*South
            + TradeLayoffs*Vetspc + Base + Rep2PartyVoteShare
            + EVANRATE + netMLS
            + median_age
            + lPopulation
            + IncomePC
            + Unemployment
            + Percent_Black
            + as.factor(Year)
            #+ State
            , data = 
              y3)
summary(lm1.1)

# Standard Errors clustered at county-level
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y3,lm1.1,y3$GeoFips)

cl4.1 <- cl[, 1:4]
cl4.1 # Table 4, Model 1 results

# 2.  Black
lm1.2 <- lm(Enrollment ~ lag1Enrollment
            + TradeLayoffs*South
            + TradeLayoffs*Percent_Black + Base + Rep2PartyVoteShare
            + EVANRATE + netMLS
            + median_age
            + lPopulation
            + IncomePC
            + Unemployment
            + Vetspc
            + as.factor(Year)
            #+ State
            , data = 
              y3)
summary(lm1.2)

# Standard Errors clustered at county-level
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y3,lm1.2,y3$GeoFips)

cl4.2 <- cl[, 1:4]
cl4.2 # Table 4, Model 2 results


# 3. Incomepc
lm1.3 <- lm(Enrollment ~ lag1Enrollment
            + TradeLayoffs*South
            + TradeLayoffs*IncomePC + Base + Rep2PartyVoteShare + Vetspc
            + EVANRATE + netMLS
            + median_age
            + lPopulation
            + Unemployment
            + Percent_Black
            + as.factor(Year)
            #+ State
            , data = 
              y3)
summary(lm1.3)

# Standard Errors clustered at county-level
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y3,lm1.3,y3$GeoFips)

cl4.3 <- cl[, 1:4]
cl4.3 # Table 4, Model 3 results


# 4. Median Age
lm1.4 <- lm(Enrollment ~ lag1Enrollment
            + TradeLayoffs*South
            + TradeLayoffs*median_age + Base + Rep2PartyVoteShare + Vetspc
            + EVANRATE + netMLS
            + lPopulation
            + IncomePC
            + Unemployment
            + Percent_Black
            + as.factor(Year)
            #+ State
            , data = 
              y3)
summary(lm1.4)

# Standard Errors clustered at county-level
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y3,lm1.4,y3$GeoFips)

cl4.4 <- cl[, 1:4]
cl4.4 # Table 4, Model 4 results

# 5. Youth
lm1.5 <- lm(Enrollment ~ lag1Enrollment
            + TradeLayoffs*South
            + TradeLayoffs*YOUTH + Base + Rep2PartyVoteShare + Vetspc
            + EVANRATE + netMLS
            + lPopulation
            + IncomePC
            + Unemployment
            + Percent_Black
            + as.factor(Year)
            #+ State
            , data = 
              y3)
summary(lm1.5)

# Standard Errors clustered at county-level
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y3,lm1.5,y3$GeoFips)

cl4.5 <- cl[, 1:4]
cl4.5 # Table 4, Model 5 results


# 6. Party ID
lm1.6 <- lm(Enrollment ~ lag1Enrollment
            + TradeLayoffs*South
            + TradeLayoffs*Rep2PartyVoteShare + Base + median_age + Vetspc
            + EVANRATE + netMLS
            + lPopulation
            + IncomePC
            + Unemployment
            + Percent_Black
            + as.factor(Year)
            #+ State
            , data = 
              y3)
summary(lm1.6)

# Standard Errors clustered at county-level
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y3,lm1.6,y3$GeoFips)

cl4.6 <- cl[, 1:4]
cl4.6 # Table 4, Model 6 results

# 7. Army Bases
lm1.7 <- lm(Enrollment ~ lag1Enrollment
            + TradeLayoffs*South
            + TradeLayoffs*Base + Rep2PartyVoteShare + median_age + Vetspc
            + EVANRATE + netMLS
            + lPopulation
            + IncomePC
            + Unemployment
            + Percent_Black
            + as.factor(Year)
            #+ State
            , data = 
              y3)
summary(lm1.7)

# Standard Errors clustered at county-level
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y3,lm1.7,y3$GeoFips)

cl4.7 <- cl[, 1:4]
cl4.7 # Table 4, Model 7 results

# Table 5: Union Density and the South

# Model 1
lm1.1 <- lm(Enrollment ~ lag1Enrollment
            + TradeLayoffs + Membership + South
            + Base + Rep2PartyVoteShare
            + EVANRATE + netMLS
            + median_age
            + lPopulation
            + IncomePC
            + Unemployment
            + Percent_Black
            + as.factor(Year)
            #+ State
            , data = 
              y3)
summary(lm1.1)

# Standard Errors clustered at county-level
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y3,lm1.1,y3$GeoFips)

cl5.1 <- cl[, 1:4]
cl5.1 # Table 5, Model 1 results

# Model 2

lm1.2 <- lm(Enrollment ~ lag1Enrollment
            + TradeLayoffs*South + Membership
            + Base + Rep2PartyVoteShare
            + EVANRATE + netMLS
            + median_age
            + lPopulation
            + IncomePC
            + Unemployment
            + Percent_Black
            + as.factor(Year)
            #+ State
            , data = 
              y3)
summary(lm1.2)

# Standard Errors clustered at county-level
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y3,lm1.2,y3$GeoFips)

cl5.2 <- cl[, 1:4]
cl5.2 # Table 5, Model 2 results

# Model 3

lm1.3 <- lm(Enrollment ~ lag1Enrollment
            + TradeLayoffs*Membership + South
            + Base + Rep2PartyVoteShare
            + EVANRATE + netMLS
            + median_age
            + lPopulation
            + IncomePC
            + Unemployment
            + Percent_Black
            + as.factor(Year)
            #+ State
            , data = 
              y3)
summary(lm1.3)

# Standard Errors clustered at county-level
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y3,lm1.3,y3$GeoFips)

cl5.3 <- cl[, 1:4]
cl5.3 # Table 5, Model 3 results

# Model 4

lm1.4 <- lm(Enrollment ~ lag1Enrollment
            + TradeLayoffs*Membership 
            + Base + Rep2PartyVoteShare
            + EVANRATE + netMLS
            + median_age
            + lPopulation
            + IncomePC
            + Unemployment
            + Percent_Black
            + as.factor(Year)
            + State
            , data = 
              y3)
summary(lm1.4)

# Standard Errors clustered at county-level
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y3,lm1.4,y3$GeoFips)

cl5.4 <- cl[, 1:4]
cl5.4 # Table 5, Model 4 results

# Model 5

lm1.5 <- lm(Enrollment ~ lag1Enrollment
            + TradeLayoffs*Membership + TradeLayoffs*South 
            + Base + Rep2PartyVoteShare
            + EVANRATE + netMLS
            + median_age
            + lPopulation
            + IncomePC
            + Unemployment
            + Percent_Black
            + as.factor(Year)
            #+ State
            , data = 
              y3)
summary(lm1.5)

# Standard Errors clustered at county-level
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y3,lm1.5,y3$GeoFips)

cl5.5 <- cl[, 1:4]
cl5.5 # Table 5, Model 5 results

# ROBUSTNESS TESTS

# Table 6: Difference-in-Differences

x <- y3[,1:59]
x <- na.omit(x)

# Model 1
lm1.5 <- lm(dEnrollment ~  
              dTradeLayoffs*South
            + dlPopulation + dnetMLS
            + dIncomePC
            + dUnemployment
            + dPercent_Black
            + as.factor(Year)
            # + State
            , data = 
              x)
summary(lm1.5)

cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(x,lm1.5,x$GeoFips)

cl6.1 <- cl[, 1:4]
cl6.1 # Table 6, Model 1 Results

# Model 2 
lm1.4 <- lm(dEnrollment ~  
              dTradeLayoffs*South
            + dTradeLayoffs*West
            + dTradeLayoffs*MidWest
            + dlPopulation + dnetMLS
            + dIncomePC
            + dUnemployment
            + dPercent_Black
            + as.factor(Year)
            # + State
            , data = 
              x)
summary(lm1.4)

# Standard Errors clustered at county-level
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(x,lm1.4,x$GeoFips)

cl6.2 <- cl[, 1:4]
cl6.2 # Table 6, Model 2 Results

# Table 7: Negative Binomial

# Model 1

m1 <- glm.nb(Enrollment ~ TradeLayoffs + netMLS 
             + EVANRATE
             + lPopulation
             + IncomePC
             + Unemployment + Base + Rep2PartyVoteShare
             + Percent_Black + Vetspc + median_age
             + as.factor(Year)
             + State, data=y3[y3$Enrollment > 0,])
summary(m1)

# robust SEs
cov.m1 <- vcovHC(m1, type="HC0")
std.err <- sqrt(diag(cov.m1))
r.est7.1 <- cbind(Estimate= coef(m1), "Robust SE" = std.err,
                "Pr(>|z|)" = 2 * pnorm(abs(coef(m1)/std.err), lower.tail=FALSE),
                LL = coef(m1) - 1.96 * std.err,
                UL = coef(m1) + 1.96 * std.err)
r.est7.1 # Table 7, Model 1 Results

# Model 2

m2 <- glm.nb(Enrollment ~ TradeLayoffs*South + netMLS 
             + EVANRATE
             + lPopulation
             + IncomePC
             + Unemployment + Base + Rep2PartyVoteShare
             + Percent_Black + Vetspc + median_age
             + as.factor(Year)
             #+ State
             , data=y3[y3$Enrollment > 0,])
summary(m2)

# robust SEs
cov.m2 <- vcovHC(m2, type="HC0")
std.err <- sqrt(diag(cov.m2))
r.est7.2 <- cbind(Estimate= coef(m2), "Robust SE" = std.err,
                "Pr(>|z|)" = 2 * pnorm(abs(coef(m2)/std.err), lower.tail=FALSE),
                LL = coef(m2) - 1.96 * std.err,
                UL = coef(m2) + 1.96 * std.err)
r.est7.2 # Table 7, Model 2 Results

# Model 3 

m3 <- glm.nb(Enrollment ~ TradeLayoffs*South + TradeLayoffs*median_age + netMLS 
             + EVANRATE
             + lPopulation
             + IncomePC
             + Unemployment + Base + Rep2PartyVoteShare
             + Percent_Black + Vetspc
             + as.factor(Year)
             #+ State
             , data=y3[y3$Enrollment > 0,])
summary(m3)

# robust SEs
cov.m3 <- vcovHC(m3, type="HC0")
std.err <- sqrt(diag(cov.m3))
r.est7.3 <- cbind(Estimate= coef(m3), "Robust SE" = std.err,
                "Pr(>|z|)" = 2 * pnorm(abs(coef(m3)/std.err), lower.tail=FALSE),
                LL = coef(m3) - 1.96 * std.err,
                UL = coef(m3) + 1.96 * std.err)
r.est7.3 #  Table 7, Model 3 Results


# Table 8: Recruitment and Placebo Test

# Model 1 - Relative Military Pay
lm1.6 <- lm(Enrollment ~ lag1Enrollment + RMCvECI
            + TradeLayoffs*South + netMLS
            + EVANRATE
            + lPopulation
            + IncomePC
            + Unemployment + Base + Rep2PartyVoteShare
            + Percent_Black + Vetspc + median_age
            #+ as.factor(Year)
            #+ State
            , data = 
              y3)
summary(lm1.6)

# Standard Errors clustered at county-level
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y3,lm1.6,y3$GeoFips)

cl8.1 <- cl[, 1:4]
cl8.1 # Table 8, Model 1 Results

# Model 2 - Reruiter Cost
y4 <- y3[,c(1:42,62)]
y4 <- na.omit(y4)

lm1.7 <- lm(Enrollment ~ lag1Enrollment + Recruiter_Cost
            + TradeLayoffs*South + netMLS
            + EVANRATE
            + lPopulation
            + IncomePC
            + Unemployment + Base + Rep2PartyVoteShare
            + Percent_Black + Vetspc + median_age
            #+ as.factor(Year)
            #+ State
            , data = 
              y4)
summary(lm1.7)

# Standard Errors clustered at county-level
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y4,lm1.7,y4$GeoFips)

cl8.2 <- cl[, 1:4]
cl8.2 # Table 8, Model 2 Results

# Model 3 - Military Deaths

y5 <- y3[,c(1:42,66)]
y5 <- na.omit(y5)

lm1.8 <- lm(Enrollment ~ lag1Enrollment + Military_Deaths
            + TradeLayoffs*South + netMLS
            + EVANRATE
            + lPopulation
            + IncomePC
            + Unemployment + Base + Rep2PartyVoteShare
            + Percent_Black + Vetspc + median_age
            #+ as.factor(Year)
            #+ State
            , data = 
              y5)
summary(lm1.8)

# Standard Errors clustered at county-level
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y5,lm1.8,y5$GeoFips)

cl8.3 <- cl[, 1:4]
cl8.3 # Table 8, Model 3 Results

# Model 4 - Placebo Test

y3 <- y3[,1:42]

# Generate 1 and 2 year leads of trade-related job losses

library(zoo)
library(sandwich)
library(plm)
y3 <- y3[order(y3$GeoFips, y3$Year),]
xpan <- pdata.frame(y3,c("GeoFips", "Year"))
y3$leadTradeLayoffs <- lead(xpan$TradeLayoffs)

y3 <- y3[order(y3$GeoFips, y3$Year),]
xpan <- pdata.frame(y3,c("GeoFips", "Year"))
y3$lead2TradeLayoffs <- lead(xpan$leadTradeLayoffs)

y4 <- na.omit(y3)

lm1.p <- lm(Enrollment ~ lag1Enrollment 
            + leadTradeLayoffs + South
            + EVANRATE
            + lPopulation
            + IncomePC
            + Unemployment + + Base + Rep2PartyVoteShare
            + Percent_Black + Vetspc + median_age
            + as.factor(Year)
            #+ State
            , data = 
              y4)
summary(lm1.p)

# Standard Errors clustered at county-level
cl   <- function(dat,fm, cluster){
  require(sandwich, quietly = TRUE)
  require(lmtest, quietly = TRUE)
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- fm$rank
  dfc <- (M/(M-1))*((N-1)/(N-K))
  uj  <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat=crossprod(uj)/N)
  coeftest(fm, vcovCL) }
require(foreign)

cl <- cl(y4,lm1.p,y4$GeoFips)

cl8.4 <- cl[, 1:4]
cl8.4 # Table 8, Model 4 Results

# Table 3 (Results)
cl3.1
cl3.2
cl3.3

# Table 4 (Results)
cl4.1
cl4.2
cl4.3
cl4.4
cl4.5
cl4.6
cl4.7

# Table 5 (Results)
cl5.1
cl5.2
cl5.3
cl5.4
cl5.5

# Table 6 (Results)
cl6.1
cl6.2

# Table 7 (Results)
r.est7.1
r.est7.2
r.est7.3

# Table 8 (Results)
cl8.1
cl8.2
cl8.3
cl8.4

